#include "CubeMap.hpp"

namespace zzz{
map<zuint,CubemapPosList*> Cubemapbase::PosListMap;

const CubemapPosList& Cubemapbase::GetPosList(zuint n)
{
  map<zuint,CubemapPosList*>::iterator muci=PosListMap.find(n);
  if (muci==PosListMap.end())
  {
    CubemapPosList *newposlist=new CubemapPosList;
    newposlist->reserve(6*n*n);
    const float allarea=4.0/n/n;
    const float halfn=n/2.0;
    for (zuint i=0; i<6; i++) for (zuint j=0; j<n; j++) for (zuint k=0; k<n; k++)
    {
      CubemapPos newpos;
      newpos.direction=CartesianDirCoordf(CubeCoord((CUBEFACE)i,k,j,n));  //should be i,k,j , so that it is correspond to ToIndex() order
                                        //the must-be-correct way is newposlist[CubeCoord((CUBEFACE)i,j,k,n).ToIndex()].direction=CartesianDirCoordf(CubeCoord((CUBEFACE)i,j,k,n))
                                        //however this must be quicker
      float rsquare=1+
        (halfn-j-0.5)*(halfn-j-0.5)/halfn/halfn+
        (halfn-k-0.5)*(halfn-k-0.5)/halfn/halfn;
      float costheta=1.0/sqrt(rsquare);
      float omiga=allarea*costheta/rsquare;
      newpos.solidangle=omiga;
      newposlist->push_back(newpos);
    }
    PosListMap[n]=newposlist;
    return *newposlist;
  }
  return *(muci->second);
}

map<Cubemapbase::SHBaseCubemapDesc,Cubemap<float>*> Cubemapbase::SHBaseCubemap;
const Cubemap<float>& Cubemapbase::GetSHBaseCubemap(int l,int m,int size)
{
  SHBaseCubemapDesc x;
  x[0]=l; x[1]=m; x[2]=size;
  map<SHBaseCubemapDesc,Cubemap<float>*>::iterator msi=SHBaseCubemap.find(x);
  if (msi==SHBaseCubemap.end())
  {
    Cubemapf *newcube=new Cubemapf(size);
    const CubemapPosList &poslist=GetPosList(size);
    int len=size*size*6;
    for (int i=0; i<len; i++)
    {
      SphericalCoordf sc(poslist[i].direction);
      newcube->v[i]=SHBase(l,m,sc.theta,sc.phi)*poslist[i].solidangle;
    }
    SHBaseCubemap[x]=newcube;
    return *newcube;
  }
  return *(msi->second);
}

const Cubemap<float>* Cubemapbase::GetSHBaseCubemapPointer(int l,int m,int size)
{
  SHBaseCubemapDesc x;
  x[0]=l; x[1]=m; x[2]=size;
  map<SHBaseCubemapDesc,Cubemap<float>*>::iterator msi=SHBaseCubemap.find(x);
  if (msi==SHBaseCubemap.end())
  {
    Cubemapf *newcube=new Cubemapf(size);
    const CubemapPosList &poslist=GetPosList(size);
    int len=size*size*6;
    for (int i=0; i<len; i++)
    {
      SphericalCoordf sc(poslist[i].direction);
      newcube->v[i]=SHBase(l,m,sc.theta,sc.phi)*poslist[i].solidangle;
    }
    SHBaseCubemap[x]=newcube;
    return newcube;
  }
  return msi->second;
}

float Cubemapbase::SHBase(int l,int m,float theta,float phi)
{
  const float sqrt2 = sqrt(2.0);
  if(m==0) return K(l,0)*P(l,m,cos(theta));
  else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
  else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
}

float Cubemapbase::K(int l,int m)
{
  float temp = ((2.0*l+1.0)*Factorial(l-m)) / (4.0*PI*Factorial(l+m));
  return sqrt(temp);
}

float Cubemapbase::P(int l,int m,float x)
{
  float pmm = 1.0;
  if(m>0) {
    float somx2 = sqrt((1.0-x)*(1.0+x));
    float fact = 1.0;
    for(int i=1; i<=m; i++) {
      pmm *= (-fact) * somx2;
      fact += 2.0;
    }
  }
  if(l==m) return pmm;
  float pmmp1 = x * (2.0*m+1.0) * pmm;
  if(l==m+1) return pmmp1;
  float pll = 0.0;
  for(int ll=m+2; ll<=l; ++ll) {
    pll = ((2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm) / (ll-m);
    pmm = pmmp1;
    pmmp1 = pll;
  }
  return pll;
}


}
